#####Violence in the Eye of the Beholder: Does Identity Shape Civilians' Perceptions towards Peacekeepers? #####

#Authors: Michael Yekple & Zlatin Mitkov

##Code for Replicating the Results##

###The ANCOVA visualization package uses an old version of GGplot2.
###Please install the correct version of ggplot by running lines 9 and 10. 
require(remotes)
install_version("ggplot2", version = "3.4.0", repos = "http://cran.us.r-project.org")
library(ggplot2)


options(scipen=999)

library(psych)
library(FSA)
library(Rmisc)
library(car)
library(multcompView)
library(lsmeans)
library(rcompanion)
library(emmeans)
library(dplyr) 
library(jtools)
library(stargazer)
library(lme4)
library(huxtable)
library(flextable)
library(dotwhisker)
library(sandwich)
library(broom.mixed)
library(ggpubr)
library(tidyverse)
library(tidymodels)
library(broom)
library(AICcmodavg)
library(multcomp)
library(haven)


#import the dataset 

UN_Data1 <- read_dta("UN_Data1.dta")


##Descriptive Analysis 

UN_Data1$familysafe_num<-as.numeric(UN_Data1$familysafe)

m1<-lm(un_responseto~ familysafe_num,data=UN_Data1)
summary(m1)

m2<-lm(fardc_responseto~ familysafe_num,data=UN_Data1)
summary(m2)


m3<-lm(un_responseto~ norm_security_index,data=UN_Data1)
summary(m3)

m4<-lm(fardc_responseto~ norm_security_index,data=UN_Data1)
summary(m4)


m5<-lm(un_responseto~ ethnic_violence_1,data=UN_Data1)
summary(m5)


m6<-lm(fardc_responseto~ ethnic_violence_1,data=UN_Data1)
summary(m6)



plot_summs(m1,m2,m3, m4, m5, m6, inner_ci_level = .9, 
           model.names = c("Approval for MONUSCO1", "Approval for FARDC2", 
                           "Approval for MONUSCO3", "Approval for FARDC4", 
                           "Approval for MONUSCO5", "Approval for FARDC6"), 
           vline = geom_vline(xintercept = 0, colour = "grey60",linetype=2)) %>%
  relabel_predictors(c(familysafe_num = "Family Insecurity",                       
                       norm_security_index = "Perception of Insecurity",
                       ethnic_violence_1 = "Experince of Ethnic Violence"))+
  theme_bw() + xlab("OLS Coefficient Estimates") + ylab("") +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  ggtitle("Security Provision") +
  theme(plot.title = element_text(face="bold"),
        legend.background = element_rect(colour="grey80"),
        legend.title = element_blank())


##Factor 

UN_Data1$treatments = factor(UN_Data1$Groups,
                             levels=unique(UN_Data1$Groups))


##Ordered factor 

UN_Data1$vignettes = factor(UN_Data1$Groups,
                            ordered = TRUE,
                            levels = c("1", "2", "3", "4", "5"))


UN_Data1$Calculategroup_2 = factor(UN_Data1$Calculategroup_1,
                                   levels=unique(UN_Data1$Calculategroup_1))



#############Main Paper ##########################




#DV1 Approval of MONUSCO Hypothesis 1  

Summarize(dv_approvalmonusco ~ treatments,
          data=UN_Data1,
          digits=3)

ancova_model <- aov(dv_approvalmonusco ~ treatments + 
                      age + education + 
                      defensiveviolence + 
                      norm_socialdominance_index + 
                      norm_security_index, data = UN_Data1)
Anova(ancova_model, type="III")

ancova_model <- aov(dv_approvalmonusco ~ treatments + age + 
                      education + defensiveviolence + 
                      norm_socialdominance_index + 
                      norm_security_index, data = UN_Data1)
postHocs <- glht(ancova_model, linfct = mcp(treatments = "Tukey"))
summary(postHocs)


##Data Visualization DV1 Approval of MONUSCO 



Sum = groupwiseMean(dv_approvalmonusco ~ vignettes,
                    data   = UN_Data1,
                    conf   = 0.95,
                    digits = 3,
                    na.rm = T,
                    traditional = F,
                    percentile  = T)


ggplot(Sum,                
       aes(x = vignettes,
           y = Mean)) +
  geom_errorbar(aes(ymin = Percentile.lower,
                    ymax = Percentile.upper),
                width = 0.05,
                size  = 0.5) +
  geom_point(shape = 15,
             size  = 4) +
  theme_bw() +
  theme(axis.title   = element_text(face  = "bold")) +
  ylab("Support for Monusco's Actions") + 
  xlab("Treatment Groups")



#Subseting for Dv1 for Hutu 1 

Summarize(dv_approvalmonusco ~ treatments,
          data=subset(UN_Data1, UN_Data1$Calculategroup_2==1),
          digits=3)

ancova_model20 <- aov(dv_approvalmonusco ~ treatments + age + 
                        education + defensiveviolence + 
                        norm_socialdominance_index + 
                        norm_security_index, data=subset(UN_Data1, 
                                                         UN_Data1$Calculategroup_2==1))
Anova(ancova_model20, type="III")

ancova_model20 <- aov(dv_approvalmonusco ~ treatments + age + education + 
                        defensiveviolence + norm_socialdominance_index + 
                        norm_security_index, data=subset(UN_Data1, 
                                                         UN_Data1$Calculategroup_2==1))
postHocs20 <- glht(ancova_model20, linfct = mcp(treatments = "Tukey"))
summary(postHocs20)




##Subseting for DV1 for Hutu 2 


Summarize(dv_approvalmonusco ~ treatments,
          data=subset(UN_Data1, UN_Data1$Calculategroup_2==2),
          digits=3)

ancova_model30 <- aov(dv_approvalmonusco ~ treatments + age + 
                        education + defensiveviolence + 
                        norm_socialdominance_index + 
                        norm_security_index, data=subset(UN_Data1, 
                                                         UN_Data1$Calculategroup_2==2))
Anova(ancova_model30, type="III")


ancova_model30 <- aov(dv_approvalmonusco ~ treatments + age + 
                        education + defensiveviolence + 
                        norm_socialdominance_index + 
                        norm_security_index, data=subset(UN_Data1, 
                                                         UN_Data1$Calculategroup_2==2))
postHocs30 <- glht(ancova_model30, linfct = mcp(treatments = "Tukey"))
summary(postHocs30)


##Subseting for DV1 for Tutsi 3 

Summarize(dv_approvalmonusco ~ treatments,
          data=subset(UN_Data1, UN_Data1$Calculategroup_2==3),
          digits=3)

ancova_model40 <- aov(dv_approvalmonusco ~ treatments + age + 
                        education + defensiveviolence + 
                        norm_socialdominance_index + 
                        norm_security_index, data=subset(UN_Data1, 
                                                         UN_Data1$Calculategroup_2==3))
Anova(ancova_model40, type="III")


ancova_model40 <- aov(dv_approvalmonusco ~ treatments + age + 
                        education + defensiveviolence + 
                        norm_socialdominance_index + 
                        norm_security_index, data=subset(UN_Data1, 
                                                         UN_Data1$Calculategroup_2==3))
postHocs40 <- glht(ancova_model40, linfct = mcp(treatments = "Tukey"))
summary(postHocs40)


##Subseting for DV1 for Nande 4 

Summarize(dv_approvalmonusco ~ treatments,
          data=subset(UN_Data1, UN_Data1$Calculategroup_2==4),
          digits=3)

ancova_model50 <- aov(dv_approvalmonusco ~ treatments + age + 
                        education + defensiveviolence + 
                        norm_socialdominance_index + 
                        norm_security_index, data=subset(UN_Data1, 
                                                         UN_Data1$Calculategroup_2==4))
Anova(ancova_model50, type="III")

ancova_model50 <- aov(dv_approvalmonusco ~ treatments + age + 
                        education + defensiveviolence + 
                        norm_socialdominance_index + 
                        norm_security_index, data=subset(UN_Data1, 
                                                         UN_Data1$Calculategroup_2==4))
postHocs50 <- glht(ancova_model50, linfct = mcp(treatments = "Tukey"))
summary(postHocs50)






##Hypothesis H2

#DV MONUSCO Approval Treatment 1 answers comparison between those who believe that group violence 
#is justified and those that do not

UN_Data1$defense_grouptreatment1 = factor(UN_Data1$defense_group1,
                                          levels=unique(UN_Data1$defense_group1))

Summarize(dv_approvalmonusco ~ defense_grouptreatment1,
          data=UN_Data1,
          digits=3)

modelOLS1<-lm(dv_approvalmonusco~defense_grouptreatment1 + age + 
                education + defensiveviolence + 
                norm_socialdominance_index + 
                norm_security_index, data = UN_Data1)
summary(modelOLS1)



##DV MONUSCO Approval Treatment 2 answers comparison between those who believe that group violence 
#is justified and those that do not 

UN_Data1$defense_grouptreatment2 = factor(UN_Data1$defense_group2,
                                          levels=unique(UN_Data1$defense_group2))

Summarize(dv_approvalmonusco ~ defense_grouptreatment2,
          data=UN_Data1,
          digits=3)

modelOLS2<-lm(dv_approvalmonusco~defense_grouptreatment2 + age + 
                education + defensiveviolence + 
                norm_socialdominance_index + 
                norm_security_index, data = UN_Data1)
summary(modelOLS2)



##DV MONUSCO Approval Treatment 3 answers comparison between those who believe that group violence 
#is justified and those that do not 

UN_Data1$defense_grouptreatment3 = factor(UN_Data1$defense_group3,
                                          levels=unique(UN_Data1$defense_group3))

Summarize(dv_approvalmonusco ~ defense_grouptreatment3,
          data=UN_Data1,
          digits=3)

modelOLS3<-lm(dv_approvalmonusco~defense_grouptreatment3 + age + 
                education + defensiveviolence + 
                norm_socialdominance_index + 
                norm_security_index, data = UN_Data1)
summary(modelOLS3)



##DV H2 MONUSCO Approval Treatment 4 answers comparison between those who believe that group violence is justified and those that do not 

UN_Data1$defense_grouptreatment4 = factor(UN_Data1$defense_group4,
                                          levels=unique(UN_Data1$defense_group4))

Summarize(dv_approvalmonusco ~ defense_grouptreatment4,
          data=UN_Data1,
          digits=3)

modelOLS4<-lm(dv_approvalmonusco~defense_grouptreatment4 + age + 
                education + defensiveviolence + 
                norm_socialdominance_index + 
                norm_security_index, data = UN_Data1)
summary(modelOLS4)



##DV H2 MONUSCO Approval Treatment 5 answers comparison between those who believe that group violence 
#is justified and those that do not 


UN_Data1$defense_grouptreatment5 = factor(UN_Data1$defense_group5,
                                          levels=unique(UN_Data1$defense_group5))

Summarize(dv_approvalmonusco ~ defense_grouptreatment5,
          data=UN_Data1,
          digits=3)

modelOLS5<-lm(dv_approvalmonusco~defense_grouptreatment5 + age + 
                education + defensiveviolence + 
                norm_socialdominance_index + 
                norm_security_index, data = UN_Data1)
summary(modelOLS5)



#Visualization of the results from Hypothesis 2  



plot_summs(modelOLS1,modelOLS2,modelOLS3, modelOLS4, modelOLS5, inner_ci_level = .9,
           omit.coefs = c("age", "education"),
           model.names = c("Approval for MONUSCO1", "Approval for FARDC2", 
                           "Approval for MONUSCO3", "Approval for FARDC4", "Approval for MONUSCO5"),
           vline = geom_vline(xintercept = 0, colour = "grey60",linetype=2)) %>%
  relabel_predictors(c(defensiveviolence = "Use of Defensive Violence",
                       norm_socialdominance_index = "Social Dominance",
                       norm_security_index = "Perception of Insecurity",
                       defense_grouptreatment10 = "T1 Group Violence is Justified",                       
                       defense_grouptreatment21 = "T2 Group Violence is Justified",
                       defense_grouptreatment31 = "T3 Group Violence is Justified",
                       defense_grouptreatment41 = "T4 Group Violence is Justified",
                       defense_grouptreatment51 = "T5 Group Violence is Justified"))+
  theme_bw() + xlab("OLS Coefficient Estimates") + ylab("") +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  ggtitle("Support for MONUSCO's Actions") +
  theme(plot.title = element_text(face="bold"),
        legend.background = element_rect(colour="grey80"),
        legend.title = element_blank())



#DV Morality Hypothesis 3 

Summarize(dv_morality ~ treatments,
          data=UN_Data1,
          digits=3)

ancova_model333 <- aov(dv_morality ~ treatments + age + 
                         education + defensiveviolence + 
                         norm_socialdominance_index + 
                         norm_security_index, data = UN_Data1)
Anova(ancova_model333, type="III")

ancova_model333 <- aov(dv_morality ~ treatments + age + 
                         education + defensiveviolence + 
                         norm_socialdominance_index + 
                         norm_security_index, data = UN_Data1)
postHocs333 <- glht(ancova_model333, linfct = mcp(treatments = "Tukey"))
summary(postHocs333)


##Data Visualization H3 

Sum = groupwiseMean(dv_morality ~ vignettes,
                    data   = UN_Data1,
                    conf   = 0.95,
                    digits = 3,
                    na.rm = T,
                    traditional = F,
                    percentile  = T)

ggplot(Sum,                
       aes(x = vignettes,
           y = Mean)) +
  geom_errorbar(aes(ymin = Percentile.lower,
                    ymax = Percentile.upper),
                width = 0.05,
                size  = 0.5) +
  geom_point(shape = 15,
             size  = 4) +
  theme_bw() +
  theme(axis.title   = element_text(face  = "bold")) +
  ylab("The Morality of the Violence") + 
  xlab("Treatment Groups")




#Subseting for Hutu 1 

Summarize(dv_morality ~ treatments,
          data=subset(UN_Data1, UN_Data1$Calculategroup_2==1),
          digits=3)

ancova_mode200 <- aov(dv_morality ~ treatments + age + 
                        education + defensiveviolence + 
                        norm_socialdominance_index + 
                        norm_security_index, data=subset(UN_Data1, 
                                                         UN_Data1$Calculategroup_2==1))
Anova(ancova_mode200, type="III")

ancova_mode200 <- aov(dv_morality ~ treatments + age + 
                        education + defensiveviolence + 
                        norm_socialdominance_index + 
                        norm_security_index, data=subset(UN_Data1, 
                                                         UN_Data1$Calculategroup_2==1))
postHocs200 <- glht(ancova_mode200, linfct = mcp(treatments = "Tukey"))
summary(postHocs200)



##Subseting for Hutu 2 

Summarize(dv_morality ~ treatments,
          data=subset(UN_Data1, UN_Data1$Calculategroup_2==2),
          digits=3)

ancova_model201 <- aov(dv_morality ~ treatments + age + 
                         education + defensiveviolence + 
                         norm_socialdominance_index + 
                         norm_security_index, data=subset(UN_Data1, 
                                                          UN_Data1$Calculategroup_2==2))
Anova(ancova_model201, type="III")

ancova_model201 <- aov(dv_morality ~ treatments + age + 
                         education + defensiveviolence + 
                         norm_socialdominance_index + 
                         norm_security_index, data=subset(UN_Data1, 
                                                          UN_Data1$Calculategroup_2==2))
postHocs201 <- glht(ancova_model201, linfct = mcp(treatments = "Tukey"))
summary(postHocs201)



##Subsseting for Tutsi 

Summarize(dv_morality ~ treatments,
          data=subset(UN_Data1, UN_Data1$Calculategroup_2==3),
          digits=3)

ancova_model202 <- aov(dv_morality ~ treatments + age + 
                         education + defensiveviolence + 
                         norm_socialdominance_index + 
                         norm_security_index, data=subset(UN_Data1, 
                                                          UN_Data1$Calculategroup_2==3))
Anova(ancova_model202, type="III")

ancova_model202 <- aov(dv_morality ~ treatments + age + education + 
                         defensiveviolence + norm_socialdominance_index + 
                         norm_security_index, data=subset(UN_Data1, 
                                                          UN_Data1$Calculategroup_2==3))
postHocs202 <- glht(ancova_model202, linfct = mcp(treatments = "Tukey"))
summary(postHocs202)




##Subseting for Nande 

Summarize(dv_morality ~ treatments,
          data=subset(UN_Data1, UN_Data1$Calculategroup_2==4),
          digits=3)

ancova_model203 <- aov(dv_morality ~ treatments + age + 
                         education + defensiveviolence + 
                         norm_socialdominance_index + 
                         norm_security_index, data=subset(UN_Data1, 
                                                          UN_Data1$Calculategroup_2==4))
Anova(ancova_model203, type="III")

ancova_model203 <- aov(dv_morality ~ treatments + age + education + 
                         defensiveviolence + norm_socialdominance_index + 
                         norm_security_index, data=subset(UN_Data1, 
                                                          UN_Data1$Calculategroup_2==4))
postHocs203 <- glht(ancova_model203, linfct = mcp(treatments = "Tukey"))
summary(postHocs203)





##Morality of Violence Treatment 1 Hypothesis 4 

#1 people think that violence to protect their group is not morally wrong  
#0 people think that violence to protect their group is morally wrong  

UN_Data1$morality_grouptreatment1 = factor(UN_Data1$morality_group1,
                                           levels=unique(UN_Data1$morality_group1))

Summarize(dv_approvalmonusco ~ morality_grouptreatment1,
          data=UN_Data1,
          digits=3)

modelOLS27<-lm(dv_approvalmonusco~morality_grouptreatment1 + age + 
                 education + defensiveviolence + 
                 norm_socialdominance_index + 
                 norm_security_index, data = UN_Data1)
summary(modelOLS27)


##DV Morality of Violence Treatment 2
#1 people think that violence to protect their group is not morally wrong  
#0 people think that violence to protect their group is morally wrong  

UN_Data1$morality_grouptreatment2 = factor(UN_Data1$morality_group2,
                                           levels=unique(UN_Data1$morality_group2))

Summarize(dv_approvalmonusco ~ morality_grouptreatment2,
          data=UN_Data1,
          digits=3)

modelOLS28<-lm(dv_approvalmonusco~morality_grouptreatment2 + age + 
                 education + defensiveviolence + 
                 norm_socialdominance_index + 
                 norm_security_index, data = UN_Data1)
summary(modelOLS28)


##DV Morality of Violence Treatment 3
#1 people think that violence to protect their group is not morally wrong  
#0 people think that violence to protect their group is morally wrong  

UN_Data1$morality_grouptreatment3 = factor(UN_Data1$morality_group3,
                                           levels=unique(UN_Data1$morality_group3))

Summarize(dv_approvalmonusco ~ morality_grouptreatment3,
          data=UN_Data1,
          digits=3)

modelOLS29<-lm(dv_approvalmonusco~morality_grouptreatment3 + 
                 age + education + defensiveviolence + 
                 norm_socialdominance_index + 
                 norm_security_index, data = UN_Data1)
summary(modelOLS29)



##DV Morality of Violence Treatment 4
#1 people think that violence to protect their group is not morally wrong  
#0 people think that violence to protect their group is morally wrong  

UN_Data1$morality_grouptreatment4 = factor(UN_Data1$morality_group4,
                                           levels=unique(UN_Data1$morality_group4))

Summarize(dv_approvalmonusco ~ morality_grouptreatment4,
          data=UN_Data1,
          digits=3)

modelOLS30<-lm(dv_approvalmonusco~morality_grouptreatment4 + 
                 age + education + defensiveviolence + 
                 norm_socialdominance_index + 
                 norm_security_index, data = UN_Data1)
summary(modelOLS30)



##DV Morality of Violence Treatment 5
#1 people think that violence to protect their group is not morally wrong  
#0 people think that violence to protect their group is morally wrong  

UN_Data1$morality_grouptreatment5 = factor(UN_Data1$morality_group5,
                                           levels=unique(UN_Data1$morality_group5))

Summarize(dv_approvalmonusco ~ morality_grouptreatment5,
          data=UN_Data1,
          digits=3)

modelOLS31<-lm(dv_approvalmonusco~morality_grouptreatment5 + age + 
                 education + defensiveviolence + 
                 norm_socialdominance_index + 
                 norm_security_index, data = UN_Data1)
summary(modelOLS31)



#Visualization for H4 


plot_summs(modelOLS27,modelOLS28,modelOLS29, modelOLS30, modelOLS31, inner_ci_level = .9,
           omit.coefs = c("age", "education"),
           model.names = c("Approval for MONUSCO1", "Approval for FARDC2", 
                           "Approval for MONUSCO3", "Approval for FARDC4", "Approval for MONUSCO5"),
           vline = geom_vline(xintercept = 0, colour = "grey60",linetype=2)) %>%
  relabel_predictors(c(defensiveviolence = "Use of Defensive Violence",
                       military_force = "Use of Military Force",
                       norm_socialdominance_index = "Social Dominance",
                       norm_security_index = "Perception of Insecurity",
                       morality_grouptreatment11 = "T1 Morality and Group Violence",                       
                       morality_grouptreatment20 = "T2 Morality and Group Violence",
                       morality_grouptreatment30 = "T3 Morality and Group Violence",
                       morality_grouptreatment40 = "T4 Morality and Group Violence",
                       morality_grouptreatment51 = "T5 Morality and Group Violence"))+
  theme_bw() + xlab("OLS Coefficient Estimates") + ylab("") +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  ggtitle("The Morality of the Violence") +
  theme(plot.title = element_text(face="bold"),
        legend.background = element_rect(colour="grey80"),
        legend.title = element_blank())











#############Appendix I ##########################

#################################################################
##Balance of the treatments


ballance1 <-lm(UN1_treatment~ age + 
                 education + 
                 defensiveviolence + 
                 norm_socialdominance_index + 
                 norm_security_index, data = UN_Data1)
summary(ballance1)

ballance2 <-lm(UN2_treatment~ age + education + 
                 defensiveviolence + 
                 norm_socialdominance_index + 
                 norm_security_index, data = UN_Data1)
summary(ballance2)

ballance3 <-lm(UN3_treatment~ age + education + 
                 defensiveviolence + 
                 norm_socialdominance_index + 
                 norm_security_index, data = UN_Data1)
summary(ballance3)

ballance4 <-lm(UN4_treatment~ age + education + 
                 defensiveviolence + 
                 norm_socialdominance_index + 
                 norm_security_index, data = UN_Data1)
summary(ballance4)

ballance5 <-lm(UN5_treatment~ age + 
                 education + defensiveviolence + 
                 norm_socialdominance_index + 
                 norm_security_index, data = UN_Data1)
summary(ballance5)

stargazer(ballance1, ballance2, ballance3, ballance4, ballance5, 
          type="html",out="combined_result_2.doc",
          title="Ballance of Treatments OLS", 
          style="apsr", intercept.bottom = T,intercept.top = F,digits=3, 
          label="tab:randcheck", dep.var.labels.include=FALSE, 
          column.labels=c("Treatment 1", "Treatment 2", "Treatment 3", 
                          "Treatment 4", "Treatment 5"), 
          covariate.labels=c("Age", "Education", "Defensiv Violence", 
                             "Social Dominance", "Perceived Security"), 
          star.cutoffs=c(0.10, 0.05,0.01))




###############################################

##Auxiliary Data Results###

##Auxiliary Hypothesis 1
##Group Self-Defense  

Summarize(dv_groupdefense ~ treatments,
          data=UN_Data1,
          digits=3)

ancova_model60 <- aov(dv_groupdefense ~ treatments + age + 
                        education + defensiveviolence + 
                        norm_socialdominance_index + norm_security_index, data = UN_Data1)
Anova(ancova_model60, type="III")

ancova_model60 <- aov(dv_groupdefense ~ treatments + age + 
                        education + 
                        defensiveviolence + 
                        norm_socialdominance_index + 
                        norm_security_index, data = UN_Data1)
postHocs60 <- glht(ancova_model60, linfct = mcp(treatments = "Tukey"))
summary(postHocs60)



#Data Visualization Group Self Defense 

Sum = groupwiseMean(dv_groupdefense ~ vignettes,
                    data   = UN_Data1,
                    conf   = 0.95,
                    digits = 3,
                    na.rm = T,
                    traditional = F,
                    percentile  = T)


ggplot(Sum,                
       aes(x = vignettes,
           y = Mean)) +
  geom_errorbar(aes(ymin = Percentile.lower,
                    ymax = Percentile.upper),
                width = 0.05,
                size  = 0.5) +
  geom_point(shape = 15,
             size  = 4) +
  theme_bw() +
  theme(axis.title   = element_text(face  = "bold")) +
  ylab("Support for Group Self-Defense") + 
  xlab("Treatment Groups")



#Subseting for Hutu 1 

Summarize(dv_groupdefense ~ treatments,
          data=subset(UN_Data1, UN_Data1$Calculategroup_2==1),
          digits=3)

ancova_model70 <- aov(dv_groupdefense ~ treatments + 
                        age + education + 
                        defensiveviolence + 
                        norm_socialdominance_index + 
                        norm_security_index, data=subset(UN_Data1, 
                                                         UN_Data1$Calculategroup_2==1))
Anova(ancova_model70, type="III")

ancova_model70 <- aov(dv_groupdefense ~ treatments + age + 
                        education + 
                        defensiveviolence + 
                        norm_socialdominance_index + 
                        norm_security_index, data=subset(UN_Data1, 
                                                         UN_Data1$Calculategroup_2==1))
postHocs70 <- glht(ancova_model70, linfct = mcp(treatments = "Tukey"))
summary(postHocs70)



##Subseting for Hutu 2 

Summarize(dv_groupdefense ~ treatments,
          data=subset(UN_Data1, UN_Data1$Calculategroup_2==2),
          digits=3)

ancova_model80 <- aov(dv_groupdefense ~ treatments + age + 
                        education + 
                        defensiveviolence + 
                        norm_socialdominance_index + 
                        norm_security_index, data=subset(UN_Data1, 
                                                         UN_Data1$Calculategroup_2==2))
Anova(ancova_model80, type="III")

ancova_model80 <- aov(dv_groupdefense ~ treatments + age + 
                        education + 
                        defensiveviolence + 
                        norm_socialdominance_index + 
                        norm_security_index, data=subset(UN_Data1, 
                                                         UN_Data1$Calculategroup_2==2))
postHocs80 <- glht(ancova_model80, linfct = mcp(treatments = "Tukey"))
summary(postHocs80)



##Subseting for Tutsi 3 

Summarize(dv_groupdefense ~ treatments,
          data=subset(UN_Data1, UN_Data1$Calculategroup_2==3),
          digits=3)

ancova_model90 <- aov(dv_groupdefense ~ treatments + age + 
                        education + 
                        defensiveviolence + 
                        norm_socialdominance_index + 
                        norm_security_index, data=subset(UN_Data1, 
                                                         UN_Data1$Calculategroup_2==3))
Anova(ancova_model90, type="III")

ancova_model90 <- aov(dv_groupdefense ~ treatments + age + 
                        education + 
                        defensiveviolence + 
                        norm_socialdominance_index + 
                        norm_security_index, data=subset(UN_Data1, 
                                                         UN_Data1$Calculategroup_2==3))
postHocs90 <- glht(ancova_model90, linfct = mcp(treatments = "Tukey"))
summary(postHocs90)



##Subseting for Nande 

Summarize(dv_groupdefense ~ treatments,
          data=subset(UN_Data1, UN_Data1$Calculategroup_2==4),
          digits=3)

ancova_model101 <- aov(dv_groupdefense ~ treatments + age + 
                         education + defensiveviolence + 
                         norm_socialdominance_index + 
                         norm_security_index, data=subset(UN_Data1, 
                                                          UN_Data1$Calculategroup_2==4))
Anova(ancova_model101, type="III")

ancova_model101 <- aov(dv_groupdefense ~ treatments + age + 
                         education + defensiveviolence + 
                         norm_socialdominance_index + 
                         norm_security_index, data=subset(UN_Data1, 
                                                          UN_Data1$Calculategroup_2==4))
postHocs101 <- glht(ancova_model101, linfct = mcp(treatments = "Tukey"))
summary(postHocs101)



#########################


##Auxiliary Hypothesis 2
##Group Violence Treatment 1 

#1 people think it that group violence is justified
#0 people do not think that group violence is justified  

UN_Data1$defense_grouptreatment1 = factor(UN_Data1$defense_group1,
                                          levels=unique(UN_Data1$defense_group1))

Summarize(dv_morality ~ defense_grouptreatment1,
          data=UN_Data1,
          digits=3)

modelOLS16<-lm(dv_morality~defense_grouptreatment1 + age + education + 
                 defensiveviolence + norm_socialdominance_index + 
                 norm_security_index, data = UN_Data1)
summary(modelOLS16)



##Group Violence Treatment 2
#1 people think it that group violence is justified
#0 people do not think that group violence is justified  


UN_Data1$defense_grouptreatment2 = factor(UN_Data1$defense_group2,
                                          levels=unique(UN_Data1$defense_group2))

Summarize(dv_morality ~ defense_grouptreatment2,
          data=UN_Data1,
          digits=3)

modelOLS17<-lm(dv_morality~defense_grouptreatment2 + age + education + 
                 defensiveviolence + norm_socialdominance_index + 
                 norm_security_index, data = UN_Data1)
summary(modelOLS17)



##Group Violence Treatment 3
#1 people think it that group violence is justified
#0 poeple do not think that group violence is justified  

UN_Data1$defense_grouptreatment3 = factor(UN_Data1$defense_group3,
                                          levels=unique(UN_Data1$defense_group3))

Summarize(dv_morality ~ defense_grouptreatment3,
          data=UN_Data1,
          digits=3)

modelOLS18<-lm(dv_morality~defense_grouptreatment3 + age + education + 
                 defensiveviolence + norm_socialdominance_index + 
                 norm_security_index, data = UN_Data1)
summary(modelOLS18)



##Group Violence Treatment 4
#1 people think it that group violence is justified
#0 people do not think that group violence is justified  

UN_Data1$defense_grouptreatment4 = factor(UN_Data1$defense_group4,
                                          levels=unique(UN_Data1$defense_group4))

Summarize(dv_morality ~ defense_grouptreatment4,
          data=UN_Data1,
          digits=3)

modelOLS19<-lm(dv_morality~defense_grouptreatment4 + age + 
                 education + defensiveviolence + 
                 norm_socialdominance_index + 
                 norm_security_index, data = UN_Data1)
summary(modelOLS19)



##Group Violence Treatment 5
#1 people think it that group violence is justified
#0 people do not think that group violence is justified  


UN_Data1$defense_grouptreatment5 = factor(UN_Data1$defense_group5,
                                          levels=unique(UN_Data1$defense_group5))

Summarize(dv_morality ~ defense_grouptreatment5,
          data=UN_Data1,
          digits=3)

modelOLS20<-lm(dv_morality~defense_grouptreatment5 + 
                 age + education + defensiveviolence + 
                 norm_socialdominance_index + 
                 norm_security_index, data = UN_Data1)
summary(modelOLS20)


#######################

##Auxiliary Hypothesis 3
##DV Excessive Violence and its morality Treatment 1

#1 people think that violent act can be attributed to the group's leadership  
#0 people think that the violent act cannot be attributed to the group's leadership  

UN_Data1$escessive_grouptreatment1 = factor(UN_Data1$escessive_group1,
                                            levels=unique(UN_Data1$escessive_group1))

Summarize(dv_morality ~ escessive_grouptreatment1,
          data=UN_Data1,
          digits=3)

modelOLS32<-lm(dv_morality~escessive_grouptreatment1 + age + 
                 education + defensiveviolence + 
                 norm_socialdominance_index + 
                 norm_security_index, data = UN_Data1)
summary(modelOLS32)


## Excessive Violence and its morality Treatment 2
#1 people think that violent act can be attributed to the group's leadership  
#0 people think that the violent act cannot be attributed to the group's leadership  

UN_Data1$escessive_grouptreatment2 = factor(UN_Data1$escessive_group2,
                                            levels=unique(UN_Data1$escessive_group2))

Summarize(dv_morality ~ escessive_grouptreatment2,
          data=UN_Data1,
          digits=3)

modelOLS33<-lm(dv_morality~escessive_grouptreatment2 + age + 
                 education + defensiveviolence + 
                 norm_socialdominance_index + 
                 norm_security_index, data = UN_Data1)
summary(modelOLS33)



## Excessive Violence and its morality Treatment 3 
#1 people think that violent act can be attributed to the group's leadership  
#0 people think that the violent act cannot be attributed to the group's leadership  

UN_Data1$escessive_grouptreatment3 = factor(UN_Data1$escessive_group3,
                                            levels=unique(UN_Data1$escessive_group3))

Summarize(dv_morality ~ escessive_grouptreatment3,
          data=UN_Data1,
          digits=3)

modelOLS34<-lm(dv_morality~escessive_grouptreatment3 + age + 
                 education + defensiveviolence + 
                 norm_socialdominance_index + 
                 norm_security_index, data = UN_Data1)
summary(modelOLS34)



## Excessive Violence and its morality Treatment 4 
#1 people think that violence to protect their group is not morally wrong  
#0 people think that violence to protect their group is morally wrong  

UN_Data1$escessive_grouptreatment4 = factor(UN_Data1$escessive_group4,
                                            levels=unique(UN_Data1$escessive_group4))

Summarize(dv_morality ~ escessive_grouptreatment4,
          data=UN_Data1,
          digits=3)

modelOLS35<-lm(dv_morality~escessive_grouptreatment4 + age + 
                 education + defensiveviolence + 
                 norm_socialdominance_index + 
                 norm_security_index, data = UN_Data1)
summary(modelOLS35)



## Excessive Violence and its morality Treatment 5 
#1 people think that violent act can be attributed to the group's leadership  
#0 people think that the violent act cannot be attributed to the group's leadership  

UN_Data1$escessive_grouptreatment5 = factor(UN_Data1$escessive_group5,
                                            levels=unique(UN_Data1$escessive_group5))

Summarize(dv_morality ~ escessive_grouptreatment5,
          data=UN_Data1,
          digits=3)

modelOLS36<-lm(dv_morality~escessive_grouptreatment5 + age + 
                 education + defensiveviolence + 
                 norm_socialdominance_index + 
                 norm_security_index, data = UN_Data1)
summary(modelOLS36)



#######################


##Auxiliary Hypothesis 4
##DV excessive Violence 

Summarize(dv_exessiveviolence ~ treatments,
          data=UN_Data1,
          digits=3)

ancova_model222 <- aov(dv_exessiveviolence ~ treatments + age + 
                         education + defensiveviolence + 
                         norm_socialdominance_index + 
                         norm_security_index, data = UN_Data1)
Anova(ancova_model222, type="III")

ancova_model222 <- aov(dv_exessiveviolence ~ treatments + age + 
                         education + defensiveviolence + 
                         norm_socialdominance_index + 
                         norm_security_index, data = UN_Data1)
postHocs222 <- glht(ancova_model222, linfct = mcp(treatments = "Tukey"))
summary(postHocs222)



##Data Visualization H4 

Sum = groupwiseMean(dv_exessiveviolence ~ vignettes,
                    data   = UN_Data1,
                    conf   = 0.95,
                    digits = 3,
                    na.rm = T,
                    traditional = F,
                    percentile  = T)

ggplot(Sum,                
       aes(x = vignettes,
           y = Mean)) +
  geom_errorbar(aes(ymin = Percentile.lower,
                    ymax = Percentile.upper),
                width = 0.05,
                size  = 0.5) +
  geom_point(shape = 15,
             size  = 4) +
  theme_bw() +
  theme(axis.title   = element_text(face  = "bold")) +
  ylab("Attributing the Violence to the group's Leaders") + 
  xlab("Treatment Groups")




#Subseting for Hutu 1 

Summarize(dv_exessiveviolence ~ treatments,
          data=subset(UN_Data1, UN_Data1$Calculategroup_2==1),
          digits=3)

ancova_mode160 <- aov(dv_exessiveviolence ~ treatments + age + 
                        education + 
                        defensiveviolence + 
                        norm_socialdominance_index + 
                        norm_security_index, data=subset(UN_Data1, UN_Data1$Calculategroup_2==1))
Anova(ancova_mode160, type="III")

ancova_mode160 <- aov(dv_exessiveviolence ~ treatments + age + 
                        education + defensiveviolence + 
                        norm_socialdominance_index + 
                        norm_security_index, data=subset(UN_Data1, UN_Data1$Calculategroup_2==1))
postHocs160 <- glht(ancova_mode160, linfct = mcp(treatments = "Tukey"))
summary(postHocs160)



##Subseting for Hutu 2 


Summarize(dv_exessiveviolence ~ treatments,
          data=subset(UN_Data1, UN_Data1$Calculategroup_2==2),
          digits=3)

ancova_model170 <- aov(dv_exessiveviolence ~ treatments + age + 
                         education + defensiveviolence + 
                         norm_socialdominance_index + 
                         norm_security_index, data=subset(UN_Data1, UN_Data1$Calculategroup_2==2))
Anova(ancova_model170, type="III")

ancova_mode170 <- aov(dv_exessiveviolence ~ treatments + age + 
                        education + defensiveviolence + 
                        norm_socialdominance_index + 
                        norm_security_index, data=subset(UN_Data1, UN_Data1$Calculategroup_2==2))
postHocs170 <- glht(ancova_mode170, linfct = mcp(treatments = "Tukey"))
summary(postHocs170)



##Subseting for Tutsi 3

Summarize(dv_exessiveviolence ~ treatments,
          data=subset(UN_Data1, UN_Data1$Calculategroup_2==3),
          digits=3)

ancova_model180 <- aov(dv_exessiveviolence ~ treatments + age +
                         education + defensiveviolence + 
                         norm_socialdominance_index + 
                         norm_security_index, data=subset(UN_Data1, UN_Data1$Calculategroup_2==3))
Anova(ancova_model180, type="III")

ancova_model180 <- aov(dv_exessiveviolence ~ treatments + age + 
                         education + defensiveviolence + 
                         norm_socialdominance_index + 
                         norm_security_index, data=subset(UN_Data1, UN_Data1$Calculategroup_2==3))
postHocs180 <- glht(ancova_model180, linfct = mcp(treatments = "Tukey"))
summary(postHocs180)



##Subseting for Nande  

Summarize(dv_exessiveviolence ~ treatments,
          data=subset(UN_Data1, UN_Data1$Calculategroup_2==4),
          digits=3)

ancova_model190 <- aov(dv_exessiveviolence ~ treatments + age + 
                         education + defensiveviolence + 
                         norm_socialdominance_index + 
                         norm_security_index, data=subset(UN_Data1, UN_Data1$Calculategroup_2==4))
Anova(ancova_model190, type="III")

ancova_model190 <- aov(dv_exessiveviolence ~ treatments + age + 
                         education + defensiveviolence + 
                         norm_socialdominance_index + 
                         norm_security_index, data=subset(UN_Data1, UN_Data1$Calculategroup_2==4))
postHocs190 <- glht(ancova_model190, linfct = mcp(treatments = "Tukey"))
summary(postHocs190)



####################

##Auxiliary Hypothesis 5
##Proportionality of the Violence


Summarize(dv_proportionality_num ~ treatments,
          data=UN_Data1,
          digits=3)

ancova_model111 <- aov(dv_proportionality ~ treatments + age + 
                         education + defensiveviolence + 
                         norm_socialdominance_index + 
                         norm_security_index, data = UN_Data1)
Anova(ancova_model111, type="III")

ancova_model111 <- aov(dv_proportionality ~ treatments + age + education + 
                         defensiveviolence + 
                         norm_socialdominance_index + 
                         norm_security_index, data = UN_Data1)
postHocs111 <- glht(ancova_model111, linfct = mcp(treatments = "Tukey"))
summary(postHocs111)



##Data Visualization

Sum = groupwiseMean(dv_proportionality ~ vignettes,
                    data   = UN_Data1,
                    conf   = 0.95,
                    digits = 3,
                    na.rm = T,
                    traditional = F,
                    percentile  = T)


ggplot(Sum,                
       aes(x = vignettes,
           y = Mean)) +
  geom_errorbar(aes(ymin = Percentile.lower,
                    ymax = Percentile.upper),
                width = 0.05,
                size  = 0.5) +
  geom_point(shape = 15,
             size  = 4) +
  theme_bw() +
  theme(axis.title   = element_text(face  = "bold")) +
  ylab("Perception of Proportionality") + 
  xlab("Treatment Groups")



#Subseting for Hutu 1 

UN_Data1$dv_proportionality_num = as.numeric(UN_Data1$dv_proportionality)

Summarize(dv_proportionality_num ~ treatments,
          data=subset(UN_Data1, UN_Data1$Calculategroup_2==1),
          digits=3)

ancova_mode120 <- aov(dv_proportionality_num ~ treatments + age + 
                        education + 
                        defensiveviolence + 
                        norm_socialdominance_index + 
                        norm_security_index, data=subset(UN_Data1, UN_Data1$Calculategroup_2==1))
Anova(ancova_mode120, type="III")

ancova_mode120 <- aov(dv_proportionality_num ~ treatments + age + 
                        education + defensiveviolence + 
                        norm_socialdominance_index + 
                        norm_security_index, data=subset(UN_Data1, UN_Data1$Calculategroup_2==1))
postHocs120 <- glht(ancova_mode120, linfct = mcp(treatments = "Tukey"))
summary(postHocs120)




##Subseting for Hutu 2 

Summarize(dv_proportionality_num ~ treatments,
          data=subset(UN_Data1, UN_Data1$Calculategroup_2==2),
          digits=3)

ancova_model130 <- aov(dv_proportionality_num ~ treatments + age + 
                         education + defensiveviolence + 
                         norm_socialdominance_index + 
                         norm_security_index, data=subset(UN_Data1, UN_Data1$Calculategroup_2==2))
Anova(ancova_model130, type="III")

ancova_mode130 <- aov(dv_proportionality_num ~ treatments + age + 
                        education + 
                        defensiveviolence + 
                        norm_socialdominance_index + 
                        norm_security_index, data=subset(UN_Data1, UN_Data1$Calculategroup_2==2))
postHocs130 <- glht(ancova_mode130, linfct = mcp(treatments = "Tukey"))
summary(postHocs130)



##Subseting for Tutsi 

Summarize(dv_proportionality_num ~ treatments,
          data=subset(UN_Data1, UN_Data1$Calculategroup_2==3),
          digits=3)

ancova_model140 <- aov(dv_proportionality_num ~ treatments + age + 
                         education + defensiveviolence + 
                         norm_socialdominance_index + 
                         norm_security_index, data=subset(UN_Data1, UN_Data1$Calculategroup_2==3))
Anova(ancova_model140, type="III")

ancova_model140 <- aov(dv_proportionality_num ~ treatments + age + 
                         education + defensiveviolence + 
                         norm_socialdominance_index + 
                         norm_security_index, data=subset(UN_Data1, UN_Data1$Calculategroup_2==3))
postHocs140 <- glht(ancova_model140, linfct = mcp(treatments = "Tukey"))
summary(postHocs140)



##Subseting for Nande 

Summarize(dv_proportionality_num ~ treatments,
          data=subset(UN_Data1, UN_Data1$Calculategroup_2==4),
          digits=3)

ancova_model150 <- aov(dv_proportionality_num ~ treatments + age +
                         education + defensiveviolence + 
                         norm_socialdominance_index + 
                         norm_security_index, data=subset(UN_Data1, UN_Data1$Calculategroup_2==4))
Anova(ancova_model150, type="III")

ancova_model150 <- aov(dv_proportionality_num ~ treatments + age + 
                         education + 
                         defensiveviolence + 
                         norm_socialdominance_index + 
                         norm_security_index, data=subset(UN_Data1, UN_Data1$Calculategroup_2==4))
postHocs150 <- glht(ancova_model150, linfct = mcp(treatments = "Tukey"))
summary(postHocs150)




####End of the Script ####








